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1.  Introduction 

The  Department  of  Defense  (DoD)  has  an  urgent 
need  to  both  understand  and  predict  the  effects  of  high- 
altitude  (z  > 10  km)  turbulence  (HAT)  on  systems  that 
operate  at  or  propagate  through  those  altitudes.  Examples 
of  systems  affected  by  HAT  include  surveillance  aircraft 
such  as  the  U-2  and  the  unmanned  Global  Hawk, 
developing  weapons  systems  such  as  the  Airborne  Laser 
(ABL),  and  developing  communication  systems  such  as 
the  Transformational  Communication  Satellites  (TSAT). 
The  surveillance  aircraft  are  primarily  affected  by 
mechanical  turbulence,  which  is  the  fluctuation  of  wind 
velocity  and  is  usually  caused  by  velocity  gradients,  flow 
over  terrain,  or  convection.  High  frequency  gravity 
waves  can  affect  the  aircraft  the  same  way  and  are  caused 
by  the  same  forcing  mechanisms  as  actual  mechanical 
turbulence  and  are  considered  as  part  of  the  HAT  forecast 
problem.  The  effects  of  mechanical  turbulence  include 
upsetting  the  autopilot,  degrading  the  quality  of  the 
surveillance  observations,  and  in  the  worse  case 
jeopardizing  the  aircraft  itself.  Laser-based  systems  such 
as  the  ABL  and  TSAT  are  affected  by  optical  turbulence, 
which  is  the  fluctuation  of  the  index  of  refraction  and  is 
caused  by  mechanical  turbulence  in  the  presence  of 
temperature  gradients.  For  directed-energy  weapon 
systems  such  as  the  ABL,  optical  turbulence  causes  the 
laser  to  wander,  spread,  and  scintillate  which  will  degrade 
its  power,  result  in  increased  dwell  time  to  destroy  a 
target,  and  decrease  the  effective  range.  For  laser-based 
communication  systems  such  as  TSAT,  optical  turbulence 
can  result  in  significant  data  dropouts. 

To  address  the  issue  of  predicting  high-altitude 
turbulence  and  its  effects  on  Air  Force  systems,  multiple 
DoD  Challenge  Projects  have  been  carried  out  over  the 
past  few  years  to  use  high  resolution  numerical 


simulations  to  better  understand  turbulence  that  occurs  in 
the  stably  stratified  atmosphere  above  the  planetary 
boundary  layer.  Specifically,  two  approaches  have  been 
used.  First,  since  it  is  generally  acknowledged  that 
synoptic  (Ax  ~ 100  km)  and  mesoscale  (1  km  < Ax  < 100 
km)  weather  patterns  are  responsible  for  generating  high- 
altitude  turbulence,  simulations  were  carried  out  nesting 
down  from  synoptic  scale  global  forecast  models  to 
mesoscale  models  to  microscale  models  (Ax  ~ O(lOOm)). 
This  allows  identification  and  better  understanding  of  the 
large  scale  trigger  mechanisms  that  cause  high  frequency 
waves  and  turbulence  at  high  altitudes.  Second,  direct 
numerical  simulations  (DNS)  were  conducted  of  limited 
domains  to  isolate  the  onset  and  evolution  of  turbulence 
for  specific  idealized  situations,  such  as  Kelvin- 
Helmholtz  (K.-H)  instabilities.  In  Section  2,  previous 
work  addressing  this  problem  will  be  summarized. 
Section  3 will  briefly  discuss  the  two  codes  that  are 
currently  being  used  to  carry  out  the  approaches  discussed 
above.  In  Section  4,  the  results  from  the  past  year  will  be 
discussed  in  detail.  In  Section  5,  the  implications  of  the 
results  presented  here  will  be  described. 

2.  Previous  Work 

Initial  DNS  of  K-H  instabilities  by  Weme  and  Fritts 
(2000)  visually  depicted  the  fact  that  optical  turbulence 
and  mechanical  turbulence  do  not  necessarily  coincide; 
mechanical  turbulence  is  large  when/where  mixing  is 
large,  while  optical  turbulence  can  be  significantly 
reduced  by  mixing.  However,  optical  turbulence  is  likely 
to  be  found  on  the  top  and  bottom  edges  of  mechanically 
mixed  regions  where  the  isothermal  layer  caused  by  the 
mixing  interacts  with  the  ambient  stably  stratified 
temperature  gradient,  thus  enhancing  the  temperature 
variances  in  those  areas.  Early  DNS  results  confirmed 
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observations  of  balloon  measurements  (Weme,  et  al., 
2001)  that  show  that  turbulence  in  the  free  atmosphere 
occurs  in  very  thin  vertical  layers  with  length  scales  on 
the  order  of  1 - 10  m.  Weme  and  Fritts  (2000)  also 
computed  2nd-order  structure  functions  for  temperature 
and  the  three  component  winds.  The  results  for  three 
wind  components  were  not  close  to  each  other,  indicating 
that  the  turbulence  was  anisotropic,  non-homogeneous 
and  thus  non-Kolmogorov.  This  was  subsequently 
confirmed  by  aircraft  observations  (Wroblewski,  et  al., 
2007). 

More  recent  simulations  Weme,  et  al.  (2005)  has 
expanded  both  the  spatial  and  temporal  domain  of  the 
K-H  instabilities.  This  work,  leveraged  with  Capability 
Application  Project  (CAP)  grants,  has  allowed  for  the 
simulation  of  the  K-H  under  different  sub-critical 
Richardson  numbers  (Ri<0.25).  This  has  greatly 
increased  the  understanding  of  the  K-H  process  in  the 
atmosphere.  Particularly  noteworthy  is  the  striking 
difference  in  the  transition  to  turbulence  under  differing 
stratifications.  However  despite  the  complexity  of  the 
transition  to  turbulence  and  its  dependence  on  Ri,  there 
are  some  simple  laws  that  help  us  in  our  parameterization 
work.  They  are  a)  the  final  layer  depth  scales  with  Ri~l/2 
and  b)  the  final  Ri  profile  tends  to  0.55. 

Work  with  the  synoptic-scale  to  microscale  nesting 
started  with  relatively  simple  idealized  cases  such  as  a jet 
stream  co-located  with  the  tropopause  (Joseph,  et  al., 
2004)  using  the  Arizona  State  University  (ASU) 
microscale  model.  Results  from  this  work  showed  that 
with  fine  spatial  resolution  in  turbulent  regimes,  one 
could  get  good  estimates  of  dissipation.  Subsequently  the 
approach  was  adapted  (Mahalov,  et  al.,  2006)  and  applied 
using  the  ASU  microscale  model  nested  within  the 
Weather  Research  and  Forecast  (WRF)  model.  This 
permitted  the  application  of  the  microscale  code  to  real 
cases.  Results  from  Mahalov,  et  al.  (2006)  showed  the 
depiction  of  mountain  waves  off  Greenland  and  the 
depiction  of  thin  layers  of  turbulence  by  the  microscale 
code  that  were  not  detected  by  the  WRF  model. 

3.  Methodology 

The  technical  approach  employed  for  the  high- 
resolution  simulations  of  atmospheric  wind  shears  uses 
the  DNS  code  named  TRIPLE  which  is  an  incompressible 
pseudo-spectral  Navier-Stokes  solver  that  uses  Fourier 
and  Chebyshev  spatial  discretizations.  It  is  a mature 
production  code  that  has  served  as  the  workhorse  for  six 
DoD  High  Performance  Computing  Modernization 
Program  (HPCMP)  and  NSF  Grand  Challenge  projects 
and  two  HPCMP  CAP  Phase  II  projects.  The  numerical 
method  advances  the  Fourier  and/or  Chebyshev 
coefficients  in  spectral  space.  Potential  truncation  errors 


are  removed  using  the  “tau”-correction  (Weme,  1995), 
and  vertical  stacking  of  Chebyshev  domains  is  possible. 
The  method  is  Nth-order  accurate  in  space,  where  N is  the 
number  of  spectral  modes.  Depending  on  the  application 
(i.e.,  wave-breaking,  shear,  convection,  etc.),  different 
spatial  representations  are  used.  Possible  boundary 
conditions  include  stress-free,  no-slip,  fixed-temperature 
or  heat-flux,  radiative,  incident  gravity  waves,  or  any 
combination  thereof.  For  shear-instability  and  reduced- 
equation  studies,  greater  computation  speeds  are  achieved 
by  employing  simpler  boundary  conditions  permitted  by 
these  problems.  For  these  cases  stress-free  top  and 
bottom  boundaries  permit  Fourier  representations  in  the 
vertical,  and  complicated  inversion  operations  can  be 
avoided. 

Time-stepping  is  conducted  using  the  mixed 
implicit/explicit  3rd-order  Runge-Kutta  (RK3)  scheme  of 
Spalart,  et  al.  (1991),  which  has  the  same  storage 
requirements  of  most  2nd-order  schemes,  allowing  us  to 
achieve  larger  numerical  domains  than  would  otherwise 
be  possible.  Nonlinear  terms  are  treated  explicitly,  while 
buoyancy  and  diffusion  are  handled  implicitly. 

Derivatives  (non-linear  terms)  are  computed  efficiently  in 
spectral  (physical)  space. 

Initial  wind-shear  solutions  (i.e.,  Weme  and  Fritts, 
2000)  were  computed  in  domains  of  size  (X,  X/3,  2X)  in 
the  (streamwise,  spanwise,  vertical)  directions,  where  X is 
the  wavelength  of  the  most  unstable  asymptotic  linear 
eigen  mode  for  the  K-H  instability.  In  the  initial  period  of 
this  Challenge  project  and  with  additional  resources  from 
the  2004  CAP  II  allocations,  expanded  runs  were 
completed  in  domains  of  size  (4L,  2X,  2X),  for  the  Ri  of 
0.05,  0.15,  and  0.20  and  Reynold’s  number  (Re)  of  2500 
(Weme  et  al.,  2005).  In  the  past  year  and  again  with  the 
help  of  additional  resources  from  the  2006  CAP  II 
program,  the  0.15  and  0.20  were  remn  with  Re  of  2,900 
and  4,000  respectively.  These  Re  were  chosen  so  that  all 
solutions  required  nearly  the  same  number  of  spectral 
modes,  which  peaked  at  (nx,  ny,  nz)  = (3,600,  1,800, 
1,800). 

For  the  meso  to  mircoscale  portion  of  this  work  the 
Advanced  Research  WRF  core  (Skamarock  et  al.,  2005) 
was  used.  The  model  uses  a time-split  integration  scheme 
where  low-frequency  modes  that  are  meteorologically 
significant  are  integrated  using  a RK3  time  integration 
scheme  and  high-frequency  acoustic  modes  are  integrated 
implicitly  in  the  vertical  over  smaller  time  steps  to 
maintain  stability.  Spatial  discretization  is  done  using  a C 
grid  staggering  where  normal  velocity  are  staggered  one- 
half  grid  length  from  the  thermodynamic  variables. 
Advection  of  vector  and  scalar  fields  is  in  the  form  of  flux 
divergence,  and  is  performed  using  the  third  order  Runge- 
Kutta  time-integration  scheme  using  a fifth-  and  third- 
order  accurate  spatial  discretization. 
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The  WRF  model  was  run  for  a case  during  the  2005 
Terrain-induced  Rotors  Experiment  (T-REX).  The 
standard  WRF  model  was  configured  with  three  nests. 
The  outermost  nest  received  its  initial  and  lateral 
boundary  conditions  from  the  European  Centre  for 
Medium-range  Weather  Forecasting’s  (ECMWF)  global 
forecast  model  which  provides  data  at  a horizontal 
resolution  of  25  km  and  91  vertical  levels  up  to  0.1  hPa. 
The  three  nested  WRF  grids  had  horizontal  resolutions  of 
15,  3,  and  1 km  for  each  succeeding  inner-nest  and  150 
vertical  levels  up  to  10  hPa.  Within  the  third  nest  a 
fourth,  microscale,  nest  is  employed.  Instead  of  using  the 
previous  Arizona  State  University  (ASU)  microscale 
code,  a new  code  has  been  developed,  based  on  the  full 
WRF  equations,  algorithms  and  subroutines.  This  WRF 
compatible  code  runs  in  a one-way  nesting  with  WRF  and 
allows  for  both  horizontal  and  vertical  nesting.  Both 
upper  vertical  and  lateral  boundary  conditions  are  relaxed 
to  the  fields  of  the  finest  third  WRF  nest,  including  the 
vertical  velocity.  Following  microscale  algorithms  from 
previous  ASU  codes,  an  adaptive  staggered  mesh  is 
implemented  in  the  vertical,  with  a denser  grid  around  the 
tropopause  and  in  the  lower  stratosphere.  The  new  WRF- 
like  code  was  run  with  a horizontal  resolution  of  333  m 
and  with  720  vertical  levels. 

4.  Results 

A.  K-H  Wind  Shear  Simulations 

Results  from  the  higher  Re  runs  show  fuller  spectra 
and  more  complete  turbulence  statistics  than  the  previous 
runs.  The  comparisons  of  the  DNS  with  aircraft 
observations  begun  last  year  (Ruggiero,  et  al.,  2006)  has 
continued  and  has  gone  beyond  the  validation  of  the  DNS 
results.  Analysis  of  high  frequency  aircraft  measurements 
by  Wroblewski,  et  al.  (2007)  has  detailed  distinct 
signatures  in  the  plots  of  temperature  and  winds  observed 
along  the  aircraft’s  flight  path.  These  plots  have  been 
compared  to  simulated  transects  through  the  K-H 
simulations.  Of  particular  interest  is  to  compare  the 
different  signatures  within  the  transects  and  identify 
particular  signatures  or  fingerprints  that  are  indicative  of  a 
particular  initial  Ri  number  and  temporal  point  in  the 
turbulence  evolution.  Figure  1 is  a comparison  of  the 
DNS  and  aircraft  data.  Both  plots  exhibit  the  classic  cliff- 
ramp  structure  in  the  potential  temperature  traces.  The 
presence  of  the  relatively  smoothed  ramps  without  steps 
indicates  that  the  aircraft  measurements  occurred  in  the 
latter  stages  of  the  cliff-ramp  evolution.  Another 
comparison  shown  in  Figure  2 shows  the  combination  of 
cliff-steps  and  cliff  ramps  that  suggest  either  an  early  time 
for  Ri=0.20  or  late  time  for  Ri=0.15.  Ultimately  the  Ri 
=0.15  scenario  was  chosen  as  the  best  match  because  of 


the  lack  of  cliffs  in  the  U wind  field.  The  real  intriguing 
point  of  these  results  is  how  given  a set  of  meteorological 
variables  (in  this  case  high-resolution  temperature  and 
wind  data),  one  can  gain  additional  information  such  as 
initial  Ri  and  age  of  the  turbulence  layer  by  comparing 
directly  with  DNS  results. 

B.  Microscale  Nesting 

The  microscale  case  study  highlighted  here  is  from 
the  TREx  Intensive  Observation  Period  (IOP)  8 that  took 
place  on  31  March- 1 April  in  the  Owens  Valley  region  of 
California.  The  primary  focus  of  the  experiment  was  to 
explore  the  structure  and  evolution  of  atmospheric  rotors 
(intense  low-level  horizontal  vortices  that  form  on  axis 
parallel  to,  and  downstream  of,  a mountain  ridge  crest)  as 
well  as  associated  phenomena  in  complex  terrain.  As 
such  it  provides  an  excellent  case  to  study  the  generation 
of  mountains  waves  and  its  associated  turbulence  at  high 
altitudes.  For  this  case  the  WRF  model  with  the  custom 
WRF  microscale  nest  was  initialized  on  3 1 March  2006  at 
00  UTC.  Figure  3 shows  a cross  section  Richardson’s 
number  and  the  square  of  the  wind  shear  across  the 
Owens  Valley  from  the  microscale  nest  valid  on  1 April 
2006  at  08  UTC.  Clearly  flow  coming  over  the  over 
mountains  from  the  west  is  causing  increased  shear  and 
turbulence  in  the  valley  itself.  But  what  is  interesting  to 
note  is  that  above  the  approximate  height  of  the  Sierra 
Nevada  mountains  to  the  west  of  the  valley,  the 
atmosphere  is  relatively  stable  up  to  10-15  km  at  which 
point  turbulence  layers  (Ri  < 0.25)  are  noted  beginning  in 
horizontal  near  the  peak  of  the  Sierra  Nevada  ridge  and 
continuing  down  wind.  Evidence  of  the  turbulence  layers 
above  10  km  on  the  microscale  nest  is  confirmed  by 
coincident  radiosonde  data  (not  shown).  The  temperature 
profile  indicates  wave  activity  occurring  above  the 
tropopause  (located  at  approximately  10.5  km  MSL) 
along  with  enhanced  wind  shear.  Estimations  of  vertical 
velocity  based  on  the  rise  rate  of  the  balloon  show  vertical 
velocities  exceeding  7 m/s  at  several  altitudes  above  the 
tropopause.  These  values  for  vertical  velocity  would  have 
a substantial  impact  on  aircraft  operating  at  these 
altitudes. 

5.  Significance  to  DoD 

The  goal  of  this  project  is  to  improve  prediction  of 
high-altitude  turbulence  and  quantify  its  effects  on  DoD 
systems.  Currently  the  microscale  nesting  work  is 
achieving  the  latter  by  providing  high  resolution 
characterization  of  the  atmosphere,  including  turbulence, 
in  support  of  testing  of  the  Airborne  Laser  that  has  begun 
this  spring.  Immediately  after  a laser  test,  the  WRF 
model  with  microscale  nesting  is  run.  The  results  from 
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the  microscale  nest  are  used  to  explicitly  calculate 
turbulence  parameters  such  as  the  Ri  and  the  refractive 
index  structure  constant,  a measure  of  optical  turbulence. 
Using  these  variables,  ABL  testers  can  characterize  the 
atmosphere  they  are  lasing  through  and  predict  future  test 
results  and  operational  capabilities. 

In  the  DNS  results  described  above  it  was 
demonstrated  that  given  a certain  set  of  atmospheric 
variables,  in  this  case  temperature  and  winds  along  a 
flight  path,  and  comparing  them  to  DNS  one  can  extract 
additional  information  about  the  atmosphere,  such  as  Ri 
and  the  age  of  the  turbulence  evolution.  This  leads  to  an 
opportunity  to  infer  the  subgrid-scale  turbulence  fields  by 
comparisons  with  DNS  results  once  a set  of 
meteorological  variables  are  obtained  from  an  operational 
mesoscale  model.  This  is  possible  because  of  progress  in 
recent  years  in  applying  Bayesian  Hierarchical  Modeling 
(BHM)  techniques  to  geophysical  applications  where  the 
character  of  BHM  is  particularly  well  suited  to  addressing 
measurement  and  model  uncertainty,  while  quantitatively 
incorporating  significant  scientific  understanding.  An 
additional  advantage  of  the  BHM  approach  is  that  it 
produces  a probabilistic  forecast  which  is  more  realistic 
for  solving  the  problem  since  mesoscale  models  cannot 
explicitly  resolve  the  turbulent  flow  and  are  therefore 
unable  to  make  a deterministic  forecast.  Currently  the 
DNS  results  along  with  the  aircraft  data  cited  above  and 
balloon  measurements  will  be  used  to  relate  conditional 
probabilities  of  the  occurrence  of  optical  turbulence  based 
on  input  from  operational  mesoscale  models  to  produce 
the  next  generation  of  optical  turbulence  prediction  for  the 
ABL. 
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Figure  1.  Transects  of  potential  temperature  and  the  u,  v,  and 
w wind  components  for  (a)  K-H  DNS  with  Ri  = 0.2  and  non- 
dimensional  time  t = 77  and  (b)  aircraft  flights  over  Wales, 
UK  6 June  2000.  For  the  DNS  plots  the  temperature  is  scaled 
by  the  product  of  vertical  potential  temperature  gradient  and 
distance  from  the  center  of  the  billow  and  the  winds  are 
scaled  by  the  mean  horizontal  wind  speed. 
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Figure  2.  Transects  of  potential  temperature  and  the  u,  v,  and 
w wind  components  for  (a)  K-H  DNS  with  Ri  = 0.15  and  non- 
dimensional  time  x = 72.6  and  (b)  aircraft  flights  north  of 
Darwin,  Australia  on  22  November  2003.  Scaling  for  the  DNS 
plots  is  the  same  as  for  Figure  1 . 
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Figure  3.  Upper  panel:  (118.56W,  117.42W)-altitude  cross- 
section  at  latitude  36.82N  for  Ri  Lower:  Cross-section  for  the 
square  of  the  shear  field  for  the  innermost  microscale 
domain  (333m  grid).  Both  plot  valid  for  8:00  UTC,  April  1, 
2006. 
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